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I. INTRODUCTION 


This thesis constitutes one part of a long range project to develop new sonar 
signal processing algorithms capable of rapidly solving sonar localization problems. At 
present, the solution of the sonar fire control problem can require a considerably 
longer time than that required for most other types of fire control problems. The long 
time required to achieve a solution can cause a significant degradation in a ship’s 
ability to avoid counterdetection, due to continuously decreasing range to the target 
during problem solution. A sonar system capable of rapid target localization without 
requiring own ship's maneuvers would greatly enhance the capabilities of our ships, 
and allow for weapon firings at longer ranges. 

The research question investigated in this thesis 1s whether or not it 1s possible 
to develop an algorithm which utilizes estimates of the local angles of arrival of 
acoustic signals incident upon a planar sonar array, knowledge of the deterministic 
effects of the ocean medium on sound propagation, and local sound-speed profiles of 
the ocean, to locate an acoustic transmitter, both in azimuthal angle and 
elevation’depression angle. In addition the model-based localization algorithm 
(hereafter referred to as the ‘localization algorithm’) was designed to provide the range, 
depth, cross range, and line-of-sight range between the acoustic transmitter and the 
receive array. 

Ray acoustics provides methods of determining ranges and propagation angles 
for transmission of acoustic signals in inhomogeneous media [Ref: 1: sect. 6.2]. The 
deterministic effects of the inhomogeneous ocean medium on acoustic signals are well 
known. From a transmitter in a Known position, it is possible to develop ray traces 
that illustrate the propagation of acoustic signals through the ocean medium. The 
intent here is to use this knowledge of sound propagation to find the transmitter’s 
location based on the estimated angles of arrival at a receive array. The estimates of 
the local angles of arrival are obtained from a frequency domain adaptive beamforming 
algorithm developed by Ziomek and Chan [Ref. 2]. This algorithm performs frequency 
domain adaptive beamforming for planar sonar arravs using a modified complex LMS 
adaptive algorithm. The algorithm generates estimates of the local angles of arrival, 


namely, the azimuthal and elevation'depression angles, of incoming acoustic signals. 


However, in a real ocean environment, these local angles of arrival do not reflect the 
true line-of-sight angles to the target. 

The localization algorithm uses the angle-of-arrival estimates, plus typical 
sound-speed profiles that are normally available to ships. In addition, it was found 
that one more piece of information 1s required to localize the target. This information 
is a model-based phase weight which ts part of a model-based signal processing 
algorithm developed by Ziomek and Blount [Ref. 3]. These phase weights are used to 
“correct for deterministic, ocean medium, phase effects due to ray bending as a signal 
propagates in the inhomogeneous ocean medium whose index of refraction (sound- 
speed profile) is a function of depth.” [Ref. 3] The phase weights were originally 
developed as part of an underwater acoustic communication problem in which receiver 
and transmitter locations were known. The form of the phase weights used will be 
presented in Chapter II. 

For the problem investigated in this thesis, transmitter location 1s unknown a 
priori and, therefore, the model-based phase weights cannot be determined in exactly 
the same manner as was done in the algorithm developed by Ziomek and Blount 
[Ref. 3]. The usefulness of the localization algorithm developed in this thesis is based 
on the availability of the model-based phase weights. The research done here is a 
feasibility study of the ability to localize an acoustic transmitter if the phase weights 
were available. The development of an algorithm to generate the model-based phase 
weights was not the subject of this research. 

The localization algorithm is limited to solving a particular class of problems. 
The localization algorithm is designed to accommodate vertical variations in sound- 
speed profile or, sound-speed profiles that are functions of depth only. Horizontal or 
range variations in sound-speed profile were not examined in this initial study because 
they constitute only a relatively small portion of ocean areas. Additionally, the 
transmitter and receiver are assumed to both be within the same sound-speed gradient. 
Finally, all case studies were conducted based on the assumption that the receiver was 
in close enough proximity to the transmitter so that the acoustic signal had not passed 
through a turning point prior to reaching the receiver. A turning point is defined as 
the point along a ray path at which the angle of propagation 1s 90 degrees with respect 
to the positive Y, or depth, axis. These three restrictions were necessary to limit the 
scope of the initial study to a size that would allow for a complete verification of the 


localization technique proposed, in the time allotted for the study. 
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Chapter II describes the theorv used to develop the localization algorithm. An 
overview of the problem and its geometry is presented, and then the computations 
leading to the algorithm are discussed. Finally, the limitations of the algorithm are 
presented. 

Chapter III consists of the computer simulation results and an explanation of the 
implementation of the theory in a computer algorithm. The output from the 
localization algorithm is compared to the known geometry, and a comparison of 
double precision versus single precision results is included. Additionally, the program 
is investigated to determine if errors develop as a function of the transmission angle 
and or depth separation. As will be shown in Chapter II, the roots of a fourth-order 
polvnomual must be determined to find the angle of transmission at the source. The 
roots for the fourth-order polynomial are found through use of an International 
Mathematical Subroutine Library (IMSL) subroutine and are verified by comparison 
With graphs of the function. These graphs also assist in determining the correct root to 
use during problem solution. 


In Chapter IV, conclusions and recommendations are presented. 


I] 


I]. THEORY 


A. PROBLEM OVERVIEW AND GEOMETRY 

Traditionally, the localization of acoustic transmitters by ships has been carried 
out by obtaining many lines of bearing to the transmitter, and comparing these with 
own ship’s motion to develop a geographic picture of the transmitter’s motion. This 
method is time consuming and usually very lacking in terms of accuracy. Due to the 
nature of the deterministic effects of the ocean medium, a great deal of information is 
contained in the angles at which acoustic energy arrives at the receiver. Extraction of 
this information from the local angles of arrival, while not a simple task in of itself, 
would greatly simplify the problem of target localization. 

As a first step in exploiting the information contained within the local angles of 
arrival, a geometry must be assumed for the problem. Figure 2.1 illustrates the general 
three-dimensional geometry used in the development of the method of target 
localization presented here. 


From Figure 2.1 the following definitions are apparent: 


* Xo» Yor Zy rectangular coordinates of the transmitter in meters. 

“Xn yer CR rectangular coordinates of the center of the receive planar 
arrav In meters. 

AX. AYsaZ cross range, depth, and Z coordinate separations, respectively, 
in meters between the transmitter and the receive array, 
where: 

AX = XR = Xp 
AY = Yr - Vo 
Ee Sa 
°e AR polar radial distance in meters from the transmitter to the 


recelve array. 
Note: AR? = AX? + AZ? 


e HDLTR polar radial distance in meters that a ray would travel in a 
homogeneous medium (constant sound-speed profile) between 
depths y, and yp based on an angle of transmission of B(¥p). 


e HIDLTX, HDLTZ distances in the X and Z directions, respectively sthateamnn 
would travel in a homogeneous medium between depths yz, 
and yp based on an angle of transmission of B(yp). 

eTiREOs line-of-sight range that a ray would travel in a homogencous 
medium between depths yg and yp based on an angle of 
transmission of B(¥p). 
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Note: HRLOS* = HDLTX2 + AY? + HDLTZ?. 


¢ RLOS line-of-sight range between the transmitter and the center of 
the receive array. 
Note: RLOS* = AX? + AY? + AZ? 


© Bly) initial angle of propagation (angle of transmission), measured 
with respect to the positive Y axis, of the acoustic signal at 
source depth y, meters. 


° Blyp) angle of arrival of incident plane wave field at depth yp 
meters. 
¢ BLOS the line-of-sight angle, as measured from the positive Y axis, 


between the transmitter and the receive array. 


Note that in Figure 2.1] the positive Y axis is defined in the direction of increasing 
depth, or in the downward direction. The coordinate system shown in Figure 2.1 is 
applicable for any relative positioning of the transmitter and receive array, even if AX, 
AY, and/or AZ are negative. Thus, the algorithm will work for any direction of arrival 
of the incident acoustic plane-wave field. 

The receive array is assumed to possess knowledge of its own depth. In addition, 
the receive arrav will have available estimates of arrival direction cosines associated 
with the local angles of arrival. These estimates are computed by the frequency 
domain adaptive beamforming algorithm. From these known quantities and 
information about the local sound-speed profile, the transmitter’s location with respect 


to the receiver shall be determined. 


B. TRANSMITTER LOCALIZATION THEORY 

Energy, whether it 1s acoustic or electromagnetic, will refract as it passes from a 
medium with index of refraction n, into a medium with index of refraction n,, provided 
that n, #n,. In this study, the ocean volume is characterized by a one-dimensional 
index of refraction (sound-speed profile) that is a function of depth. Snell’s law is 


given by (Ref. 1: p. 218}, 


sin B(v) 7 sin B(¥,) 


2) 
c(y) C(¥q) i 


where c(v) is the speed of sound in meters per second at a depth y. From Snell's law a 


ray parameter may be defined as 
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sin B(y,) _ sinBlyp)  sinB(ypp)_ 1 


b= ——_— = FO 22) 
C(¥) CYp) c(yyp) C(Yyp) 
where: 
eb is the rav parameter. 
Pp is the depth of a turning point. A turning point is defined as the point 


along a ray path at which the angle of propagation, B(y-yp), 1s equal to 
90 degrees. 


At this point Biv p) 1s Known, since the direction cosine 
VR) = cos B(yR) (2.3) 


is calculated by the frequency domain adaptive beamforming algorithm. The speed of 
sound at depth yp, denoted c(yp), is normally known aboard ship as a result of 
measurements made by onboard sonar systems. 

It is assumed that the sound-speed profile is a linear function of depth with 
constant gradient. In most areas of the ocean this 1s a good approximation if both the 
transmitter and the recelve array are in the same portion of the sound-speed profile. A 
tvpical sound-speed profile is shown in Figure 2.2. The parameter g is the constant 
gradient of the sound-speed profile in seconds™!. From the surface to about 100 meters 
a positive gradient is typically observed with a gradient g © +0.016 sec’! 
{Ref. 4: p. 30], {[Ref. 5: p. 401]. Below 100 meters a negative gradient 1s present, and in 
this example g * -0.02956 sec"!. Finally, at depths between 700 to 1500 meters 
[Ref. 4: p. 32] the gradient reverts to a positive value of g * +0.017 sec! 
{Ref 5: p. 401]. The value of g in the negative portion of the gradient was computed 
by assuming the speed of sound to be 1500 meters per second at the ocean surface and 
1475 meters per second at a depth of 1000 meters [{Ref. 6: p. 3]. A depth of 1000 
meters was chosen as the starting point of the second positive gradient. The negative 
gradient was then calculated to fit between the positive gradients. Based on the 
assumption that both the transmitter and the receive array are in the same gradient of 


the sound-speed profile, the speed of sound at depth y can be found from 


CSc BY - Yq) . (2.4) 


Speed of Sound 
0 —_—_—————— COO eee 
Range 


g = 0.016/sec Sound-Speed Profile 


100 m 


g = -0.0296/sec 
1000 m 


Depth 


g = 0.017/sec 


Figure 2.2. Typical Sound-Speed Profile. 


The radius of curvature that describes the arc of the circle followed by an 
acoustic field propagating through this mediuin is then [Ref. I: p. 237] 


Op eee! 2.5) 
© [g sin B(¥9)I g sin B(yp)I | 


All the terms on the far righthand side of Equation 2.5 are known. Equation 2.4 may 
be rewritten as expressed bi Ziomen |(Iel slo eos) 


i c(y) i C(¥q) . (2.6) 
g 


‘Theretore: 
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l l 
AY = yp + ¥p =~ lp) = ol¥9) (2.7) 


and, from equations 2.1 and 2.2, 


sin B(vR) 
c(¥p) = ——*- (2.8) 
b 
and 
sin B(yg) 
oo (2.9) 
b 
Combining equations 2.8 and 2.9 with equation 2.7 it is readily observed that, 
AY — Byp) . in B(yQ) (2.10) 
= Vp -V, = — Sin B(vp) - — sin Bly 
; R “0 bg =i R bg Yo 
Or, 
AY = yp - Yq = asin B(yp) - a sin BlyQ) (211) 
where 
| (2.12) 
a=—., vte 
bg 
The only unknowns now in equation 2.11 are Bly) and AY. Also note that 
R= |a]| = radius of curvature. (2.13) 


c 


The radial distance AR shown in Figure 2.1 can be found by utilizing the 


following equation [Ref. 1: p. 238}: 


i 


e C(Yq) ae 
ea ed (%) [cos B(yy) - cos B(y)] (2.14) 


which is the Z coordinate of a ray propagating in the YZ plane. In this thesis a more 
general class of problem is assumed so that the coordinate axes can remain fixed 
relative to the platform on which the planar array is mounted. Therefore, in a three- 


dimensional system, z and Z, are replaced by the polar coordinates r and ry to give 


tidy a AN) [cos B(v,) - cos B(y)] , CANS) 
g sinB (yq) se 


and, as a result, 


l 
AR = tp -Tp = [cos B(y,) - cos Blyp)] (2.16) 


OT 


AR = rp - fy = acos B(yq) - a cos B(yp) - (Qi 


The only unknowns in equation 2.17 are B(y,) and AR. Also, note in Figure 2.1 that if 


AX = 0, (2aS}) 
then 
AR = AZ. (219) 


At this point ray acoustics cannot provide any further information to develop a 
solution to the problem. However, a model-based phase weight for a planar sonar 
array, similar to that shown in Figure 2.3 , can be used to localize the transmitter. As 
derived by Ziomek and Blount [Ref. 3] 


® (f) = -2nfyndy + Myp(Gn) nn = -(N-1)/2,...,0,...,(N-1)/2 (2.20) 





Figure 2.3 Receive Planar Array Geometry. 


where 





, -Vpf 
CYR) 
Va = cos B(y7), (222) 
and: 
. ® (f) is the phase weight in the Y direction associated with element (m,n) 
in the receive dinave 
of is the frequency of the transmitted electrical signal. 


¢@My;p(fin) is the model-based phase weight which is related to the deterministic 
angle modulation performed by the ocean medium on the transmitted 
electrical signal as a function of depth [Ref. 3]. 


YT is the depth of the transmit array. 
dv is the interelement spacing in the Y direction associated with the 
receive array. 

Equation 2.20 describes the phase weights in the Y direction that a planar sonar 
array using the three-dimensional FFT beamformer presented by Ziomek and Blount 
[Ref. 3] would use to receive an acoustic signal transmitted from a depth yy and 
received at a depth yp. This equation can be seen to consist of two parts. The first 
portion is the term -2mfyndy which is the phase weight used in traditional beam 
steering. The second part, My,p(f,n), is further described by Ziomek and Blount 
[Ref. 3] as 


Pyplfn) = - (k(yp)/2vpiie(yy)/glin plyp + ndy)- 1) + AY} (2.23) 
where the wave number in radians per meter as a function of depth yy is 


K(y7) = 2nficlyy) , (2.24) 


f=f +k, 0k = Keke (2.25) 
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(yy) 


Nnn(Vp + ndv) = —— , P26 
and, 

aK: = Spa - Ve = ndy- : (ee 
where: 


ef. is the carrier frequency in hertz, 


ef, 1s the fundamental frequency in hertz of the finite Fourier series representation 


of the complex envelope of the transmitted electrical signal, and 


eK is the highest harmonic used in the finite Fourier series. 


The term NDH defines an index of refraction which is corrected for the distance 


that the (m,n) element in the receive array is offset in the Y direction from the center 


of the array. This compensation is provided by AY, , which computes the depth 
separation between the center of the transmut array and the element (m,n). 

When using these model-based phase weights it should be noted that yy 1s 
equivalent to yp) of Figure 2.1. Additionally, 


em Coo Py) = COS Pa) (2.28) 


Dividing equation 2.24 by 2vp yields 


K(y-+) ent l Th 
ee (228) 
2VB cy7) 2Vp CYT )¥p 
The term n’p(yp + ndy) - 1 may also be rewritten as 
. c(V ) ow 
n'p(vp + ndy)- 1 = ——L — -1 (2.30) 


c(yy) + gAY, 


and, as a result, 


ZA 


c(y-7) - ely) - gAY, 


n’ + ndv)- 1 = 21 
pvp + ndy) ce (2.31) 
eee ee e4¥ (2.32) 
n ndy) - | = —————_*—_-.. 
DVR 2 El ap) ae BANC 
Therefore, substituting equations 2.32 and 2.29 into equation 2.23 yields 
-mf — [-c(y7)AY, + c(y~)AY, + gAY? 
Do o(fn) = A Levy, + oD AY, + BAY, aii 
CY7)VpB [c(yy7) + gAY] 
Tf gAY ? 


“clypvg fey) + gAY,] 


Expanding the denominator of the second term on the right side of equation 2.34 


results in 


CY) BAY SG eQ a) ae C0 pia cama (2.35) 


From equation 2.4 it can be seen that 


clyy) + &Y¥p- yz + ndy) = clyp + ndy). (2.36) 
Therefore, 
e(v) + eA = cpanel) & (2e3eh) 


Substituting equation 2.37 and equation 2.22 into equation 2.34 gives 


-TEZAY 7 


D5 p(L0) = Se 
MD c(yr)c(yp + ndy) cos B(y7) 


(2.38) 


From equation 2.2, with yp = yo 


pp 


sin B(y-7) 


(vp) = (2.39) 
b 
and, as a result, 
Dy pH (6n) pieaNe (2.40) 
; Sa se : 
ee C(V¥p tndy) sin Bly) cos B(y-y) 
From equation 2.12, 
(2a) 
a=— ; 
bg 
or 
| 
a 
Using equation 2.41 in equation 2.40 yields 
io (fn) -TfAY (2.42) 
1) = ————— 
Ne ac(Vp + ndy) sin B(yyp) cos Bly) 
where 
AY = ae - ey = ndy = AY + ndy : (2.43) 


Mette center element of the recelve atiay is chosen.as the element at which the 


phase weight D,,p(f,n) is calculated, then n = 0, and 
ie oe RT: (2.44) 


0 


Therefore, atn = 0 


z5 


-mfAY~ 


By ptf) = ‘ac(yR) sin BI) 608 BY) (2.45) 
Now squaring both sides of equation 2.10 will result in 
Le = a*sin*B(y,) - 2a*sin Biyp) sin B(yy) + a’sin*B(yp) (2.46) 
which can be used in equation 2.45 to replace AY”. Now let 
x = sin B(y,) (2.47) 
and 
y = cos Bly). (2.48) 


The x and y defined in equations 2.47 and 2.48 are not the x and y coordinates 
related to Figure 2.1. Rather, this x and y are merely dummy variables to be used in 
the solution of equation 2.45. Due to the definitions of equations 2.47 and 2.48 a 
relationship between the variables x and y is apparent, that is, 
and, therefore, 


y= + (1- x22, (2.50) 


Next, replace yy with y, in equation 2.45, substitute equation 2.46 into equation 


2.45, and multiply equation 2.45 by ac(Vp)xy. This results in 
ac(Vp)Pyp(f0)xy = -nf [a?x? - 2a’sin B(yp)x + a’sin’B(yp)] . (2.51) 


Divide both sides of equation 2.51 by mfa’ to get 
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C(Vp) 
tfa 





PyppH(F0)xy = -x? + 2 sin Blyp)x - sin*B(yp). (eo?) 


Rewriting equation 2.52 yields 








es ER @ pp (COs - 2 sin B(yp)x + sin’B(yp) = 0 (2.53) 
Or 

Ax? + Bxy + Cx + D=0 (2.54) 
where 

ees || (Oe (2.55) 

B= “Ro, (r0), (2.56) 

nfa 

C = -2 sin B(yp), (2.57) 
and 

D = sin? B(yp). (2.58) 
Substituting equation 2.50 into equation 2.54 vields 

Ax? + Bx(1-x*)!’2 + Cx + D = 0 (2.59) 
OT, 

Ax? + Cx + D= + Bx(1- x2)? (2.60) 
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Squaring both sides of equation 2.60 gives 

(Ax* + Cx +D)* = B*x2(1 - x7) = Bex? - Btx?. (2.61) 
Expanding equation 2.61 yields 

A?xt + 2ACx? + (2AD + C?)x? + 2CDx + D? = Bx? - B2x4 (2.62) 
Or 


(A* + B2)x4 + 2ACx? + (2AD - B? +C*)x* + 2CDx + D? = 0. (2.63) 


To find the unknown, x, the roots of equation 2.63 must be computed. These 
roots will also be the roots of equation 2.59. In the computer algorithm written to 


implement this theory, the value 
F(x,y) = Ax? + Bx(1 - x2)/2 + Cx + D (2.64) 


was also calculated to verify the validity of the roots found for equation 2.63. 


Recalling that equation 2.47 defined 
x = sin B(¥p) (2.47) 
and equation 2.48 defined 
y = cos B(¥p) , (2.48) 


we see that once x and y are known they may be substituted into equations 2.11 and 
2.17 to solve for AY and AR (since Biv), the radius of curvature (a), the receive array 
depth, and the local sound-speed profile are all known). 

At this point AY, AR and B(¥9) are known. The next values to be found are AX, 
AZ, RLOS, and BLOS. Using the definitions of the direction cosines as presented by 
Ziomek-(hetariesp2.20) 


v = cos P(y), (2.65) 


26 


b= cos Cy), (2256) 


and 
W = cos ¥(vV) (267) 
W ene: 
° a(v) is the angle at a depth y measured from the positive X axis to the 
vector of interest. 
e y(y) is the angle at a depth y measured from the positive Z axis to the 


MeEClOUlO! iterest. 
Referring to Figure 2.1, the direction cosine v(y) at the transmitter depth can be 


Written as 





is) Biv) = (2.68) 
VUV = Cos ae ee _—_— P 
“0 a HRLOS 
and, as a result, 
AY 
HRLOS = . (2.69) 
VW(¥q) 
Also from Figure 2.1 it can be observed that 
HRLOS* = HDLTR? + AY?. (2.70) 


In ray acoustics, as presented by Ziomek [Ref. 1: p.223], the propagation vector 


is defined as 


ie tkyy 1 k5Z (2201) 
where 
ky = ku, ee 


2] 


ky = kyv, (2.73) 





7 = kw , (2.74) 
and 
Live 
ky = | (2.75) 
C(¥p) 


Therefore at the transmltcer 





amt 
alae UGE). (2.76) 
ol c(y-p) T 


and, at the receive array, 





KXR = u(yR). (2.77) 


2K 
c(V R) 


Additionally, for an inhomogeneous medium which has a sound-speed profile that is a 


function of depth only, it is known that [Ref. 1: p.223] 
kyR = kyr = CONStant. (2.78) 


Therefore, from equations 2.76 and 2.77, 





anf nf 579 
yp) WYR) = ae u(y) (2.79) 
so that 
uvy) = SR uyp , (2.80) 
(vy) 


Oi, 


c(Yy) 
‘t) = ——~ulyp). 2.81 
u(y 7) typ) U(YR) (2.81) 


In equation 2.81 u(yp) is supplied by the beamformer, c(yp) is known by own 
ship, and since equation 2.11 has been solved for AY, it 1s possible to use equation 2.4 
to calculate c(y7). Therefore, u(yy) becomes a known quantity. An alternate method 
of determining c(y7) would be by the use of Snell’s law, or equation 2.1, since Biyp), 
B(yy) and C(Vp) are all Known. 

Similarly [Ref. 1: p. 233], 








k7R = kzy (2.82) 
and, as a result, 
eae) = a W(¥Rp). (2.83) 
Referring to Figure 2.1 and utilizing equation 2.81 it can be seen that 
iy = COontt.,) = Ay) ae (2.84) 


ULV se, 
yp) *  HRLOS 


Therefore, since u(y,) is known from substituting yg for yp in equation 2.81, the value 
of HDLTX is given by 


HDLTX = u(y,)HRLOS . (2.85) 


Now that u(y,) and v(v)) are known from equation 2.84 and equation 2.68, it is 


possible to find W(¥Q) by use of the fact that [Ref. 1: p. 224] 


w*(y ES Ci : ¥*(¥q) : (2.86) 


9) 


pay) 


Again utilizing Figure 2.1 and the fact that 


W(¥o) = 608 1(¥q) (2.87) 


we see that 


ee HDLTZ (2.88) 
WIV a : 
"0" HRLOS 
Therefore, 
HDLTZ = w(y)) HRLOS . (2.89) 


Figure 2.4 shows the geometry of Figure 2.1 as seen by looking down into the 
XZ plane from above the transmitter’s depth. From Figure 2.4 the relationships 
between AZ, AX, and AR may be derived. 


The angle 6 in Figure 2.4 can be found from 


HDLTX 
tan 6 = = (2.90) 
HDLTZ 
Therefore, 
6 = tan !}(HDLTX/HDLTZ) . (3) 


Substituting equations 2.85 and 2.89 into equation 2.91 results in 


§ = tan”! {{u(y,)HRLOS}/[w(y,)HRLOS}} (2.92) 
so that 
6 = tan” !{u(y9)/w(yo)] - (2.93) 
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Figure 2.4 Topview of Geometry. 
For an inhomogeneous medium with a sound-speed profile that is a function of 


depth only, it can be shown that [Ref. I: p. 232] 


u(y) 


w(y) 





= constant. (2.94) 


Therefore, 
6 = tan"[u(yp) WR) (2.95) 


where u(Yp,) and W(yp) are available from the frequency domain adaptive beamformer. 
From Figure 2.4, AZ and AX are given by 


AZ = AR cos 6 (2.960) 


and 


Dal 


AX = AR cos 6. (22975 


Referring once again to Figure 2.1, RLOS 1s given by 


RLOS = (AX? + AY? + AZ?)!/? (2.98) 
Or, since 

AR? = AX? + AZ’, (2.99) 

RLOS = (AR? + AY?)!/2, (2.100) 


Finally, BLOS can be determined by using 


BLOS = cos (AY/RLOS) . (2.101) 


The equations presented in this section comprise the theory used to develop the 
model-based localization algorithm. By the use of ray acoustics and the assumption 
that the model-based phase weight is known, a closed form solution is possible for the 
localization problem. Obviously, the solution’s accuracy depends on a ship’s ability to 
correctly measure the sound-speed profile and the effects of any other local sonar 
conditions, such as shallow depths and the presence of biologics. However, in the open 
ocean, when the transmitter and receiver are located in the same gradient of the sound- 
speed profile, a reasonably accurate solution is possible. There are some limitations 
involved with the use of ray acoustics and the model-based phase weights. These 


limitations will be discussed in the next section. 


C. LIMITATIONS OF RAY ACOUSTICS SOLUTION 
1. Turning Points 
A turning point is that position along a ray path propagating through an 
inhomogeneous medium at which the angle of propagation measured with respect to 
the positive Y axis, B(y), is equal to 90 degrees. At this point the origination of the ray 


path becomes ambiguous to a receiver using the localization technique described in this 
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thesis, because there is no way of knowing how many turning points the acoustic 
signal has passed through. The turning point will cause a transmitter that is below 
(above) the receiver to appear to be above (below) the receiver. Figure 2.5 illustrates 
these two possibilities. In the case of receiver one in Figure 2.5, a turning point has 
occurred between the transmitter’s location and the receiver's location. The theory 
presented in this section would result in a calculated line-of-sight sinular to that shown 
in Figure 2.5. The acoustic signal passes through two turning points prior to reaching 
receiver two, and the resulting line-of-sight calculation would indicate that the 


transnutter 1s at a depth below receiver two. 










speed of Sound 
0 
Range 
g = 0.016/sec Sound-Speed Profile 
100 m 
g = -0.0296/sec 
1000 m @ Transmitte Calculated LOS 
fe 
2 
. Receiver One ete 
g = 0.01 7/sec a 


Figure 2.5 Turning Point Ambiguity. 


The turning point ambiguity problem is not necessarily very restrictive, 
depending on local sonar conditions. Table | Hsts the location of turning points in 
terms of AY and AR between the transmitter and receive array. The values in Table 1 


were calculated by assuming the values for B(¥y), c(¥9), and g shown in Table I, and 


GQ 
Coe 


then using equations 2.11 and 2.17 with Biyp) = 90°. These results show that for 
most angles of transmission the floor of the ocean would be reached prior to the signal 
reaching a turning point. Even at angles of transmission greater than 60 degrees the 


ranges to a turning point are quite large. 


TABLED 


DEPTH AND RANGE TO TURNING POINTS FOR A POSITIVE 
GRADIENT 


Biv) AY (km) AR (km) 
10° 412.913 492.043 
ZO 166.935 238.343 
30° 86.765 150.276 
40° 48.241] 103.424 
50° Boe | 86.765 
60° 13.449 50.063 
70° eo) Shoo 
80° 1.336 ese 
85° 0.331 1.092 


C(¥p) = 1475 m/sec g = 0.017 sec"! 





If the transmitter and receive array are located in the negative gradient 
portion of the sound-speed profile as shown in Figure 2.5, the situation becomes much 
more restrictive. Here the transmitter must transmit in the upward direction to reach a 
turning point, as opposed to the downward transmission assumed in Table 1. Table 2 
contains the results of calculations for the turning points in this region. In this case, 
the angles were only varied from 91 degrees to 100.8 degrees in order to place the 
turning point within the negative portion of the sound-speed profile of Figure 2.5. 
Even with the higher magnitude gradient used in Table 2, ranges of several thousand 
meters are achievable prior to the turning point. Note that all distances in Table 2 are 


in meters, whereas those listed in Table 1 are in kilometers. 
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TABIEE 2 


DEPTH AND RANGE TO TURNING POINTS FOR A NEGATIVE 


Biv) 
91.0° 
93.0° 
95.0° 
97.0° 
99.0° 
100,0° 
100.2° 
100.4° 
100.6° 
100.8° 


c(y,) = 1475 misec 


GRADIENT 


AY (m) 
-7.601 
-68.478 

- 190.604 
-374.729 
-621.991 
-769.765 
-801.279 
-833.451 
-§66.283 
-899.777 


g = -0.02956 sec”! 


2. Changes in Sound-Speed Profile 


AR (m) 

870.982 
2615.070 
4365.554 
6126.767 
7903.148 
8798.454 
8978.159 
9158.09] 
9338.253 
9518.650 


The transmitter and receiver must be in the same gradient of the sound-speed 


profile for the theory presented in this thesis to work. If the transmitter and receiver 


were located in different gradients of the sound-speed profile, a false location would be 


indicated due to the change in local angle of arrival. This situation is illustrated in 


Figure 2.6. 


3. Walidity of Model-Based Phase Weights 


The development of the model-based phase weights is based in part on the 


assumption presented by Ziomek [Ref. I: p.253] that if 


In-(y) - W/v*(yp) < < 1, 


then 


ky(y) = ky + ky*[n’(y) - I] (2ky) 


(2.102) 


(2.103) 


Speed of Sound 








0 
Range 
Sound-Speed Profile 
Transmitter 
Kes 
7 
Q 
MY 
Q 
Calculated LOS 
Receiver 
Figure 2.6 Changing Sound-Speed Gradient. 
wiiere 
vq) 
n(y) = zai (2.104) 
C(Yp) 
and 
anf 
Ke KY) VG aa V(¥p) - (2.105) 
(Vy) 


For some cases, such as B(y)) approaching 90 degrees, v(¥y) becomes very 
small, resulting i the criteria of equation 2.102 being violated. In these instances the 
model-based phase weights can no longer be considered valid. Computations were 


performed prior to running the test cases presented in this thesis to ensure that test 
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cases which violate equation 2.102 were identified and not misrepresented as valid test 
cases. 

In addition, the WKB approximation, which is the basis for the development 
of the model-based phase weights, becomes invalid as ky(y) approaches zero 
Pxctwlp 2ls|. Inisis the Case at a tuming point. 

4. Depth Separation of Zero Meters 

If AY = 0.0, meters the angle of transmission, Bly), and the local angle of 
arrival, B(vp), must both be equal to 90 degrees to permit the receive array to receive 
anv signal without that signal having to pass through a turning point. The algorithm 
fails here due to its invalidity at turning points and, as can be observed in equation 
2.17, because AR would always be computed as zero. Obviously, a AY = 0.0 meters 
does not necessarily imply that AR = 0.0 meters, since this condition is normally 


known as a collision. 


Ba 


Il. COMPUTER IMPLEMENTATION OF LOCALIZATION THEORY 


A. PROGRAM DESCRIPTION 

The implementation of the theory described in Chapter II was performed by 
writing the FORTRAN computer program LOCATE. LOCATE 1s designed to operate 
as a subroutine in the frequency domain adaptive beamforming algorithm developed by 
Ziomek and Chan [Ref. 2]. LOCATE contains one subroutine, PLOTER, which 
creates plots of the function described by equation 2.64. The description of LOCATE 
that follows demonstrates the relationship between the equations of Chapter II and the 
flow diagrams, however, the actual FORTRAWN statements are not presented. After 
LOCATE 1s explained, there 1s a short discussion of PLOTER. Section B discusses the 
method by which the algorithm was validated. Section C provides the actual results as 
compared to known geometries, and gives a comparison of double precision versus 
single precision results. 

1. Program LOCATE 

The program LOCATE uses as inputs the estimated direction cosines for local 
angles of arrival, model-based phase weights, and knowledge of the local sound-speed 
profile to determine AZ, cross-range (AX), depth separation (AY), and the line-of-sight 
range (RLOS) to the transmitter. Also, elevation/depression angle and azimuthal angle 
to the transmitter are provided by LOCATE. 

The elevation/depression angle, as shown in Figure 3.1], 1s defined as the 
minimum angle between the receive planar sonar array’s XZ plane and the line-of-sight 
between the transmitter and the receive array. The elevation/depression angle 1s 
defined to be positive (elevation) if the transmitter’s depth is less than the receiver's 
depth. If the transmitter is at a greater depth than the receive array the 
elevation/depression angle is negative (depression). Therefore the elevation/depression 
angle ranges in value from -90 degrees to +90 degrees. 

The azimuthal angle, as shown in Figure 3.2, is defined as the minimum angle 
between the receive planar sonar arrav’s Z axis and the line-of-sight between the 
transmitter and receive array, in the receive array’s XZ plane. The azimuthal angle 
then ranges from +180 degrees to 0 degrees for positive AX and from 0 degrees to 
-180 degrees for negative AX. | 

The inputs to the program LOCATE are defined as follows: 
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XZ plane 
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Figure 3.1 Elevation;Depresston Angle. 
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estimates of direction cosines (YP), WR), and W(Vp ), 
respectively, as calculated by the frequency domain 
adaptive beamformer. 


model-based phase weights. 
carrier frequency of the received electrical signal. 


fundamental frequency of the finite [ourter — serics 
FOPIcsSemiationme stew complexn envelope of the received 
electrical srgnal. 


gradient of local sound-speed profile. 

speed of sound at receive array depth yp. 

total number of recerve elements along the receive array’s Y 
AXIS. 

parameters used to deternune which harmonic ts to be used 
mn current calculations. 


parameter used to determine which element’s phase weight 
to use. 


All the inputs are currently available from the frequency domain adaptive 


beamforming algorithm described by Ziomek and Chan [Ref. 2], with the exception of 


Piiieereiires 3-3 and 3.4 illustrates the flow of the program LOCATE. 


Be 





negative 
azimuthal! 


angle ——— 


Receive array's 
Z axis 


positive —— 


azimuthal 
angle 


Figure 3.2. Azimuthal Angle. 


The beamforming algorithm is written in single precision FORTRAN. 
llowever, the programy LOCATE niust operate in double precision to enable ita 
develop accurate roots for equation 2.63. Therefore, the values passed to LOCATE 
from the adaptive beamforming algorithm must be converted to double precision, 
either in LOCATE, or before they are sent te ZO@ATE. In this thesis; all vatmes 
passed to LOCATE were double precision values. [For testing purposes, only the 
portions of the adaptive beamforming program which develop values required by 
LOCATE were used, along with a program entitled SOUNDRAY, which generates the 
true problem geometry. The reasons for the use of double precision and the support 
programms used in testing LOCATE are further described in Section [II.B.1. 

Once the program LOCATE is entered=gea  locp eiiparamicies 
QTEMP = 1, QTOTAL 1s established. From QUEMIP, an index © for the harniome 
of interest is chosen. This value Q is then used to determine the frequency F that will 
be used for further computations. 

To calculate the ray parameter SMB the local angle of arrival, B(yz), is first 


found by the arc cosine of VYR. Then SMB is calculated by equation 2.2, using CYR 


40 


START > 
.20P 
FOURTH ORDER POLYNOMIAL 


CALL ZRPOLY TO SOLVE FOR 
FOURTH ORDER POLYNOMIAL ROOTS 


NO [COMPUTE F(X,Y) FOR 
NEGATIVE VYR 
YES 


COMPUTE F(X,Y) FOR 
POSTIVE TE 


inemre 3.3 ProgrumgeOCATE I-lovwchart. 


4] 


RISLCT = G'VYR 
<His.cr > XS pf CHOOSE SMALLEST 
ROOT 


CHOOSE 
ROOT 


CALCULATE DELTAY 

AND DELTAR 
CALCULATE RLOS 
CALCULATE AZIMUTHAL ANGLE 
DELTA 

CALCULATE DELTAX 

AND DELTAZ 
CALCULATE ELEVATIONDEPRESSION 
ANGLE (ELEVDEP 
NO 

Q=QTOTAL? RETURN TO LOOP 


YES 


ee SUBROUTINE PLOTER 


a — 
RETURN TO LOOP Q = QTOTAL 
YES 


. 


: 


Figure 3.4 Program LOCATE Flowchart. 


and B(yp). The value of G is passed to LOCATE by the adaptive beamformer, so 
SMA, the radius of curvature, is now found by using equation 2.12. At this point, 
equations 2.55 through 2.58 are utilized to determine the coefficients A, B, C, and D. 
These coefficients are in turn used to find the coefficients of equation 2.63, which are 
stored in an array called COEFF. 

To determine the roots of equation 2.63, the double precision IMSL 
subroutine ZRPOLY is called, using the array COEFF as the input. ZRPOLY returns 
complex roots for equation 2.63 in an array called LAMBDA. In all the test cases that 
were run, the four roots in LAMBDA always consisted of two real roots and two 
complex roots. As a check of the validity of the roots, the value of F(x,y) from 
equation 2.64 was calculated. A graph of the function F(x,y), such as that shown in 
Figure 3.5, was used to determine whether to use +(1 - x2)h/2 or -(I - ree in this 
computation of F(x,y). 

The graphs indicated that for positive values of VYR the real roots are 
associated with the +(1l - coy term, while the complex roots are associated with the 
-(] - — term. This can be seen in Figure 3.5 where the curve associated with 
-(1 - x)? does not cross the F(x,y) = 0 line. The graph in Figure 3.5 only shows a 
small portion of the X axis. Test runs demonstrated that F(x,y) increases as x varies 
from the x value corresponding to the minimum value of F(x,y), in both the positive 
and negative X directions over the range 0 S x S 1. Therefore, the graphs were 
expanded in the region close to the minimum of F(x,y) to provide better resolution. 

To continue with the calculations, one of the four roots must be selected as 
the value x of equation 2.47. No logic in the theory section, however, provides any 
basis for a decision as to which root 1s correct. The complex roots were disregarded 
because they cannot equate to x in equation 2.47. In order to determine a relationship 
which would allow programming logic to select the correct root from the two real roots 
found by ZRPOLY, numerous test cases with known transmitter and receive array 
locations Were run using the four possible geometries allowed by the constraints listed 
on page six of this thesis. These geometries are: 

l. transmitter above receive array, 0° = B(y,) < 90°, G > 0. 
2. transmitter below receive array, 90° < Bly) = 180°,G > 0. 
Boum Gy <4), 


3. transmitter above receive array, 0° S B(y9) < 
4. transmitter below receive array, 90° < B(y,) S 180°,G < 0. 
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Figures 3.5 through 3.8 correspond to each of the four types of geometries 
listed above. Analysis of these graphs, along with Knowledge of the true geometry for 
each case, determined that if the product G* VYR, designated RTSLCT (root 
selection) in Figure 3.4, is negative, the largest real root is the correct value of x. If 
RTISLCT is positive, then the smallest root is the correct root. All test cases which 
Were subsequently run using this root selection logic resulted in the correct localization 
of the target. 

The root selected corresponds to x in equation 2.47 and is next used to 
calculate DELTAY (AY) and DELTAR (AR), using equations 2.11 and 2.17, 
respectively. From DELTAY and DELTAR , RLOS is computed from equation 2.100. 
The azimuthal angle is calculated next by equation 2.95, since UYR and WYR are 
known from the adaptive beamforming algorithm. Now DELTAZ (AZ) and DELTAX 
(AX) may be computed from equations 2.96 and 2.97, respectively. 

The elevation/depression angle is the last value to be computed. This is done 
by using equation 2.101, which provides BLOS. The angle BLOS is then converted to 


the elevation/depression angle by equation 3.1. 


ELEVDEP = 90) 300s (3.1) 


This elevation/depression angle is more useful than BLOS to personnel 
onboard ship because it provides a target location that is referenced to own ship’s 
horizontal plane. Note that computing ELEVDEP in this manner results in a negative 
value if BLOS > 90°, which indicates that the transmitter is below the receive array, 
and a positive value when BLOS < 90°, which implies that the transmitter is above the 
receive arrav. 

Program LOCATE next calls the subroutine PLOTER, if desired, to generate 
a plot of F(x,y) such as that shown in Figure 3.5. Once the graphing subroutine 1s 
completed, LOCATE checks the index @ togdetermine 11 the required) numbersor 
harmonics have been evaluated, and proceeds to process another harmonic if this has 
not been done. Otherwise, the program returns to the adaptive beamforming program. 

2. Subprogram PLOTER 

The purpose of subprogram PLOTER is to provide a graphic representation 
of the roots which the IMSL subroutine ZRPOLY calculates. The inputs to this 
subprogram are: | 


e A, B, C, D coefficients for equation 2.64. 
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Figure 3.6 F(x,y) for Geometry 2. 
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Figure 3.8 F(x,y) for Geometry 4. 
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gradient of the local sound-speed profile. 


cross-range, depth, and Z _— coordinate 
Separations calculated by the program 
POG IE. 

the line-of-sight angle as calculated by 


BOGE 


The values of G, DELTAX, DELTAY, and DELTAZ are printed out on the graph as 


G, AX, AY, and AZ, respectively, to provide a means of identifying the geometry of 


the case corresponding to each graph. 


The values A, B, C, D, and G are converted to single precision values prior to 
being passed from LOCATE to PLOTER, because PLOTER was written using 


DISSPLA which operates only in single precision. Due to the single precision accuracy 


of DISSPLA, plots made by PLOTER are not accurate enough to determine the roots 


of equation 2.64. However the plots do show approximately where the roots occur. 


Figure 3.9 illustrates the flow of the subroutine PLOTER. 











DETERMINE X COORDINATE FOR 
MINIMUM F(X,Y) = XMIN 


COMPUTE F(X,Y) FOR 
XMIN - 0.025 < X < XMIN + 0.025 


PLOT F(X,Y) VERSUS X 





Figure 3.9 Subprogram PLOTER Flowchart. 
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The subroutine PLOTER first computes the minimum value of F(x,y) in the 
interval 0 < x < 1 by incrementing x by 0.1 units. This minimum, XMIN, is then 
used as the center of the plot, with XMIN - 0.025 and XMIN + 0.025 as the lower 
and upper bounds of the graph. If XMIN + 0.025 2 1.0 the plot is centered at 0.975 
to avoid having the computer attempt to calculate the square root of a negative value 
of 1-x? in equation 2.64. After the plot is completed, PLOTER returns to the 
program LOCATE. 


B. ALGORITHM VALIDATION 
1. Generation of Received Signals 

The inputs listed in Section A of this chapter for the program LOCATE were 
generated through the use of two programs. The first program is titled SOUNDRAY 
and was written by Professor L. J. Ziomek at the U. S. Naval Postgraduate School, 
Monterey, California, in 1987. The second program is the subroutine PHSWGT 
developed by Ziomek and Blount [Ref. 7]. SOUNDRAY utilizes ray acoustics and 
geometry to develop feasible geometries for calculations of local angles of arrival of 
acoustic signals. The inputs to SOUNDRAY are the X, Y, and Z coordinates of the 
transmitter, the X and Y coordinates of the receive array, the initial angle of 
propagation, P(yy;), and information describing the local sound-speed profile. 
SOUNDRAY then uses equation 2.1 to determine B(yp) and equation 2.15 to calculate 
AR. From this point, geometry alone allows calculation of the RLOS and BLOS, from 
equations 2.98 and 2.101, respectively, and 


AZ = (RLOS? - AX? - AY?)!/2, (3.2) 


In addition, SOUNDRAY calculates the inputs for the subroutine PHSWGT 
and the estimates (in this case exact values) of direction cosines for the acoustic signal 
arriving at the receive array. SOUNDRAY determines the exact problem geometry, 
independent of the model-based phase weights, thereby providing the standard by 
which to judge the solutions generated by the program LOCATE. 

2. Test Case Results 
a. Double Precision LOCATE versus True Geometry 
As stated previously, there are four basic geometries that the program 
LOCATE is designed to handle. These four geometries may be summarized as: 
|. FAY, 0S Ce ae Once) 
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-AY, 90° < B(y,) S 180°, G > 0. 
+AY,0° S Bly.) < 90°,G < 0. 
4. -AY, 90° < Bly.) S 180°, G < 0. 


Other variations on these geometries are possible by using -AX and -AZ, 


tod 


but, because the sound-speed profile is assumed to be a function of depth only, the 
plane-wave field will propagate in a plane which is normal to the XZ plane [Ref. 1: p. 
234]. The result is that variations using -AX and -AZ merely change the sign of the 
solutions and not the magnitude. LOCATE was written to accommodate -AX and 
-AZ. However, for this discussion, it is sufficient to deal with +AX and +AZ and 
realize that only the sign of the answer is different when negative quantities are used. 

Tables 3, 4, 5, and 6 represent results from the four geometries mentioned 
above. The sound-speed profile of Figure 2.2 was used in these computations. The 
value of AY for each table was maintained constant and this necessitated the altering 
of AX depending on the angle B(y9) mised, ob BCyp) was close to 0 degrees or 
180 degrees, a smaller AX was required than for angles near 90 degrees. This is due to 
the fact that at angles near 0 degrees or 180 degrees, the plane-wave field reaches depth 
Yp in a much shorter AR than when BP(yp) is near 90 degrees. Since from equation 
2.99 


AR* = AX? + AZ?, (2.99) 


AX had to be kept sufficiently small to maintain AZ > 0, because we are working with 
cases of positive AX and AZ. 

As can be seen in Tables 3 through 6, the program LOCATE provides 
excellent results. The slight errors that are present are due mainly to roundoff error 
occurring in the root finding subroutine ZRPOLY. Note that the constraints 
concerning turning points have all been observed in these results. The maximum error 
for any range calculated by LOCATE in these cases was 0.345 meters. The angles 
calculated by LOCATE are not presented in tabular form because they were all 
accurate to four significant digits when compared to the true solutions. 

Some of the results in Tables 3 through 6 appear to be exact. This is not 
actually the case because the values in these tables were all rounded to the third 
decimal place. In no instance were the results of LOCATE exactly equal to the true 
solution, however, in many instances, the difference was in the fourth or fifth decimal 


place. 
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los 


253 
S05 
Som 
40° 
45° 


ai 
60° 
65° 


1) 
80° 
85° 


T = true solution 
+0.017 sec! 


b. Errors as a Function of Angle of Transmission and/or Depth Separation 
Figure 3.10 shows the error in RLOS as the depth 
separation between the transmitter and the receive array increases, with B(¥q) constant. 


There does not seem to be any relation between the error and the depth separation. 


AX (m) 


1 
50.000 
50.000 

100.000 
100.000 
100.000 
100.000 
100.000 
300.000 
300.000 
300.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 


L 
SeeOsr 
50.017 

100.033 
100.017 
100.015 
100.017 
100.010 
300.037 
300.020 
300.029 
500.131 
500.121 
500.049 
500.064 
500.038 
500.019 


TABLE 3 
LOCATE VERSUS TRUE GEOMETRY (GEOMETRY 1) 


AY (m) 


T 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 
300.000 


i 
300.000 
299 oe 
299.999 
SE 
300.000 
300.000 
Boo e 
300.000 
300.000 
300.000 
300.001 
300.001 
29289 10 
300.001 
300.001 
300.000 


L = LOCATE calculation 


(1) Depth Separation. 


The error appears to be mainly caused by roundoff. 


error in RLOS as the angle of transmission changes for four different depth 


separations. Again, it is readily observed that the depth separation has little effect on 


the size of the error. 


SV 


(2) Transmission Angle andjor Depth Separation. 


AZ (m) 


T 
17.434 
63.096 
44.287 
28.212 
141.879 
185.308 
231.794 

24.554 
alee 
308.993 
153.761 
414.580 
670.746 


1035.436 
1741.051 
5226.883 


L 
17.449 
63.118 
44.303 
93:229 
141.900 
185.340 
231.819 

24.559 
197.206 
307,0Zs 
153.801 
414.680 
670.811 


1035.569 
1741.185 
5227-050 





Figure 3.11 shows the 


B(v9) 


25° 
100° 
105° 
110° 
eS” 
120" 
IS" 
[50° 
133° 
140° 
145° 
150° 
[oom 
160° 
165° 
a: 


AX (m) 


re 

500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
400.000 
300.000 
200.000 
200.000 
200.000 
100.000 
100.000 
100.000 

50.000 

50.000 


iE 

500.060 
500.056 
500.013 
500.037 
500.052 
500.062 
400.027 
300.129 
200.065 
200.022 
200.051 
100.072 
100.062 
100.066 

50.047 

50.063 


TABLE 4 
EOCATE ERSUS TREE GEOVIETRY (GEOMETRY 2) 


AY (m) 


qT 
- 300.000 
- 300.000 
- 300.000 
-300.000 
- 300.000 
- 300.000 
-300.000 
- 300.000 
-300.000 
- 300.000 
- 300.000 
- 300.000 
- 300.000 
- 300.000 
- 300.000 
- 300.000 


E 
-299.993 
-299.985 
-299.976 
-300.000 
- 300.000 
-300.000 
- 300.000 
- 300.000 
- 300.000 
- 300.000 
-300.000 
- 300.000 
-300.000 
- 300.000 
- 300.000 
- 300.000 


AZ (m) 


ak 


2907.664 
1536.780 


971.804 
640.734 
3959290 
128.008 
147.308 
[Fioe3 
D22NGS 
151.643 
62.335 
140.799 
Dio 
43.151 
62.552 
16.779 


Ee 


2908.009 
[556.752 


Silas | 
640.782 
Boo 32 
128.024 
147.318 
LEGS 7 
111.211 
151.660 
627353 
140.901 
ios 
43.182 
Oe. 
16.804 


T = true solution L = LOCATE calculation 


+0.017 sec7! 





The error does increase as the angle B(y,) is increased above about 
60 degrees. This increase can be attributed to the behavior of the sine and cosine 
functions. Figure 3.12 shows how the sine and cosine functions behave between 0 and 
90 degrees. Above about 60 degrees, the slope of the sine function is less than 
0.01 degrees”! so that small changes in the sine cause large differences in the angle B(y). 
Also, in this region the magnitude of the slope of the cosine function is near its 


maximum. Small changes in the angle B(y) create large differences in the cosine. 
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B(yp) 


10° 
1S) 
20° 
20° 
30° 
cM 
40° 
45° 
50° 
a5" 
60° 
65° 
70° 
ee 
80° 
So" 


AX (m) 


I 

50.000 
100.000 
100.000 
100.000 
100.000 
300.000 
300.000 
300.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 


L 
49.995 
99.994 
oes 
100.000 

BoE 
300.004 
300.001 
300.009 
500.008 
500.006 
A99 979 
499.976 
B99 997 
499.995 
500.002 
500.004 


TABUES 
LOCATE VERSUS TRUE GEOMETRY (GEOMETRY 3) 


AY (im) 


If 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 


Ie 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
499.999 
499.999 
500.000 
500.001 
500.000 
500.000 
499.999 
499.998 


AZ (m) 


T 
nO? 
88.099 

150.837 
209.070 
268.783 
175.431 
288.24] 
393.838 
310.993 
494.933 
686.650 
916.323 
1221.188 
1671.182 
2426.104 
3902.854 


L 
207 
88.094 

150.830 
209.067 
20577 7 
175.434 
288.243 
393.850 
310.999 
494.940 
686.620 
Slee? 
1221S 
1671.169 
2426.112 
3902.89 1 


T = true solution L = LOCATE calculation 


G = -0.02956 sec”! 





To find AY, equation 2.11 uses the roots of equation 2.63 as 
determined by ZRPOLY. These roots correspond to sin B(yQ). The root contains 
some small errors due to roundoff which is borne out by the fact that the values of AY 
in Tables 3 through 6 contain errors on the order of 10°? meters. To find AR by using 
equation 2.17, the arc sine of the root must first be calculated. This amplifies any error 
in the root, especially when the angle 1s greater then 60 degrees as discussed previously. 
Next, the cosine of the arc sine of the root is computed, which further amplifies the 


CLOT. 
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B(yq) 


100° 
105 
110° 
Se 
120° 
PS” 
130° 
135° 
140° 
P95), 
150° 
I 
160° 
165° 
170° 


AX (m) 


T 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
500.000 
300.000 
300.000 
100.000 
100.000 
100.000 
100.000 

50.000 


iB 
499.998 
500.000 
499.983 
500.000 
499.986 
499.961 
A979 GD 
499.980 
299.990 
29919 
O77) 
oo 
99.990 
99.974 
49.978 


Paine ee) 
LOCATE VERSUS TRUE GEOMETRY (GEOMETRY 4) 


AY (m) 


T 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 


E 
- 500.000 
-500.000 
- 500.000 
- 500.000 
-500.000 
- 500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 
-500.000 


AZ (m) 


T 


3539.164 
1966.455 
1347.650 


983.983 
728.904 
D25292 
337.478 

71.388 
298.440 
185.560 
272.886 
212.229 
35508 

90.286 

Tae2N0 


L 


353 9Eos 
1966.417 
1347.606 


9330992 
728.884 
3297752 
337.455 

71339 
298.432 
185.548 
DEES 
212.201 
153.288 

90.264 

densi 


T = true solution L = LOCATE calculation 


G = -0.02956 sec’! 





Therefore, above about 60 degrees, we see these increased errors 
manifest themselves in the AR, AX, AZ, and RLOS calculations. Still, the errors seen 
in Figure 3.11 and in Tables 3 through 6 are insignificant when compared with the 
ranges in question. The angles are still accurate to four significant digits, and 
consequently, the range errors remain small. 

c. Double Precision Versus Single Precision Results 
It was found that the single precision version of ZRPOLY was not accurate 
enough to calculate the correct answers. The reason for this can be seen in Table 7 
which contains some single precision results for comparison to double precision results. 


ZRPOLY calculates the roots shown in the two right hand columns of Table 7. Even 
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Figure 3.10 Error in RLOS as a Function of Depth Separation. 
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Error in RLOS as a function of Transmission Angle and/or Depth Separation. 
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Figure 3.12 Sine and Cosine for 0 to 90 Degrees. 
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DP 
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60.41° 
65.44° 
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TABLE 7 
DOUBLE VAECISION VEKRSUs SINGLE PRECISION RESULTS 


DP 


RLOS (m) 


603.147 
Sou 
888.832 
1188.473 
1836.237 
9297.24 


SP 
RLOS (m) 
9.264 
48.291 
28.699 
30.835 
57.808 
350.602 


DP 

Root 
0.8660 
0.9063 
0.9397 
0.9659 
Oot 
0.9961 


SP 

Root 
0.8689 
0.9092 
0.9428 
0.9693 
0.9881 
0.9997 


+0.017 sec”! 





though the roots appear accurate to the second significant digit in the single precision 
results, when dealing with sines and cosines, an error in the third significant digit can 
create a fairly large error in calculating the angle B(y,)). Also, these roots are multiplied 
by the radius of curvature, a, in equation 2.11. This radius of curvature is on the order 
of 10° meters, so small errors in the roots will create large errors in the ranges 
calculated. The single precision results in Table 7 are so poor that they seem to have 
no relation to the actual answer. The double precision results for RLOS 1n Table 7 are 


accurate to within 0.1 meters of the true solution. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


The goal of this thesis was to determine if an underwater acoustic transmitter can 
be localized using ray acoustics, model-based phase weights, estimates of the local 
angles of arrival, and knowledge of the local sound-speed profile. As demonstrated in 
Chapter III, this goal is achievable and to a high degree of accuracy depending on the 
accuracy of the inputs to LOCATE. There are restrictions on the use of this 
procedure. It appears that the restrictions do not impose severe limitations on the use 
of the algorithm, and in some cases it may be possible to overcome them altogether. 

All the restrictions basiclly result in a limitation on the effective range of the 
algorithm. Even though acoustic signals may not reach their initial turning points for 
theoretical ranges in the tens or even hundreds of kilometers, the ocean is only about 
11.5 kilometers deep at its greatest depth. Therefore, the ranges shown in Table | are 
not realizable in some cases because the signal will reach the ocean floor in less range 
than it would take to reach the turning point. Additionally, underwater acoustic 
transmitters are usually limited in the depth to which they may be deployed, so that the 
angles of transmission that are associated with the greatest ranges will pass well below 
the receive array at any significant range. Still, the algorithm appears to be quite 
useable in ranges of less than 10 kilometers. This would be of a great advantage in the 
case of a transmitter whose signal is of low power, resulting in a short detection range. 
In fact, the need for an algorithm of this sort is most critical when the transmitter 1s at 
short range and its exact location and direction of motion must be resolved rapidly. 

In some instances, the limitations due to turning points may not be of much 
concern. For example, the algorithm might be used for an array located on the ocean 
floor. In this case, much longer ranges would be achievable, provided that the 
transmitter is in the same portion of the sound-speed profile as the receive array. The 
algorithm might also be of use in active sonar systems to provide more accurate range 
and depth information than 1s currently available. 

Implementation of the algorithm must include a very accurate root finding 
technique as has been discussed. Due to the sensitivity of the problem in regard to the 
sine and cosine functions, the roots need to be accurate to at least three significant 
figures. It was found that this is only possible through use of a double precision root 


finding subroutine. This, of course, causes the program to run more slowly but, 
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because the remainder of the program can still be written in single precision, it is not a 
great hinderance. 
In the future some areas requiring more study are: 


¢ Develop a method for obtaining the model-based phase weights from the 
received signals. At present, phase weights are computed based on received 
signals, however, the phase weights in the Y direction need to be separated into 
traditional phase weights and model-based phase weights. 


e Determine a method to account for the acoustic signal passing through a 
turning point prior to reaching the receive array. This would greatly extend the 
range capability of the algorithm. 


e Develop methods to identify signals that are transmitted from portions of the 
sound-speed profile other than the gradient in which the receive array 1s located. 


e Investigate the practical applications of the algorithm in varying acoustic 
conditions, particularly in regions such as near the Gulf Stream where the 
sound-speed profile is a function of depth and range. 
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